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ABSTRACT 

Using hydrodynamic simulations, we investigate the physical properties of gaseous substructures in 
barred galaxies and their relationships with the bar strength. The gaseous medium is assumed to be 
isothermal and unmagnetized. The bar potential is modeled as a Ferrers prolate with index n. To 
explore situations with differing bar strength, we vary the bar mass /bar relative to the spheroidal 
component as well as its aspect ratio 1Z. We derive expressions as functions of /bar and 1Z for the 
bar strength Qb and the radius r(Qb) where the maximum bar torque occurs. When applied to 
observations, these expressions suggest that bars in real galaxies are most likely to have /bar ~ 0.25- 
0.50 and n £ 1. Dust lanes approximately follow one of a^-orbits and tend to be more straight under 
a stronger and more elongated bar, but are insensitive to the presence of self-gravity. A nuclear ring 
of a conventional X2 type forms only when the bar is not so massive or elongated. The radius of 
an a^-type ring is generally smaller than the inner Lindblad resonance, decreases systematically with 
increasing Qb, and slightly larger when self-gravity is included. This evidences that the ring position 
is not determined by the resonance but by the amount of angular momentum loss at dust-lane shocks. 
Nuclear spirals exist only when the ring is of the X2-type and sufficiently large in size. Unlike the 
other features, nuclear spirals are transient in that they start out as being tightly-wound and weak, 
and then due to the nonlinear effect unwind and become stronger until turning into shocks, with an 
unwinding rate higher for larger Qb- The mass inflow rate to the galaxy center is found to be less 
than 0.01 M Q yr _1 for models with Qb ^ 0.2, while becoming larger than 0.1 M yr _1 when Qb £ 0.2 
and self-gravity is included. 

Subject headings: hydrodynamics — galaxies: ISM — galaxies: kinematics and dynamics — galaxies: 
nuclei — galaxies: spiral — galaxies: starburst — galaxies: structure — ISM: 
kinematics and dynamics — shock waves 
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1. introduction 

The inner parts of barred galaxies contain interesting 
gaseous substructures that in clude dust lanes, nuclear 
rings, and nuclear spi rals (e.g., |Sanders k Huntley||1976 
Sakamoto et al.||199 9; Knap en" et al.||2002| |Martini et al. 
2003a|bl |Sheth et al.|]2005HPrieto et ; al.||20051 |Martinez- 
Valpuesta et al.||2006| IComeron et al.||2009[ |20101|Maz- 
zuca et al.||2011| IHsieh et al.l 1201 lh . These substruc- 



tures are thought to form as a result of gas redistribu- 
tion i n galaxies initiated by non-axisymmetric bar torque 



|Combes k Gerin||1985||Buta|[l9 86; Shlosma n et al 
1990| IGarcia-Barreto et al.||1991| |Buta k Combes||1996| 
Combes||2001[ ). Understanding tneir formation and evo- 
lution is therefore of crucial importance to understand 
how the interstellar gas in barred galaxies is driven in- 
ward to affect star formation activities in the nuclear 
regions and potentially fueling of active galactic nuclei 
( AGN) . Since their spatial locations and shapes are likely 
to be determined by the ba r strength as well as the under 
lying rotation curve (e.g., Athanassoula||1992b 
k Martini||2006l IComeron et al.||2009| 



[1992b 
201Up, they may 



Peoples 



also be used to probe the mass distributions in barred 
galaxies that are not directly observable. 

While the bar substructures have been ob served for 
a long time (e.g., Pease 1917 Sandage 1961), it is dur 



ing recent years that high-quality data on their physical 
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properties have bee n compiled and compared with bar 
characteristics (e.g.,|Knapen et al.|2002[|Peeples k Mar-| 
tini 2006; Mazzuca et al. 20081 IComeron et al. 12009, 2010; 



lazzuca et al 



2011 



There is increasing observational 
evidence that the shape and size of the bar substruc- 
tures are det ermined primarily by the bar strength Q b . 



For instance, Knapen et al. (2002) measured the curva- 



ture angles Aa of dust lanes tor a sample of 9 barred 
spiral galaxies, and found that stronger bars favor more 
straight dust lanes, co nfirming a theoretical prediction of 
Athanassoula ( 1992b[) . By extending the sample size to 



55 galaxies, Comcr on et al. (2009) confirmed the results 
of Knapen et al. 12002), although they also found that a 
large scatter in the fit of the Aa-Qb relationship can be 
reduced if the bar ellipticity is also c onsidered in the fit. 
More recently, Comeron et al. (2010) measured the sizes 
and shapes of nuclear rings in a sample of 107 galaxies, 
finding that 'stronger bars host smaller rin gs' and that 
the ring e llipticity is in the range of 0-0.4. |Mazzuca et| 



al. (2011 1 reported from the analysis of 13 nuclear rings 
tnat the ring size is well correlated with the compactness 
of the galaxy mass distribution such that higher mass 
concentration implies a smaller ring. On the other hand, 
by analyzin g dust morphology in th e nuclear regions of 
75 galaxies, Peeples & Martini (2006) found that tightly- 
would nuclear spirals are preferentially found in galaxies 
with a weak bar. 

On the theoretical side, the formation and evolution 
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Kim et al. 



of bar substructure has been extensively studied using 
numerical simulations on hydrodynamical models (e.g 



Sanders fc Huntley||1976l |Roberts et al.|[l979l |van Al- 
bada fc Roberts||1981| |A thanassoula 1992b Pine r et"aT 
1995| Englmaicr & Gerhard 199 7||Patsis fc Ath anassoula 
2000| ]Maciejewski et al.||2002[ |Maciejewski||2004a|b| |Re- 
gan fc Teuben|20031 120041 |Ann fc Thakur|2005[|Lin et al.| 



2008 ; jThakur et al.||2009 ) and magnetized models (e.g 



at 



Kulesza-Zydzik et al. 2009 2010 Kulpa-Dybclct al. 2011 
Kim fc Stone|2012* ). In particu lar, |At hanassoula (1992b 
conhrmed the early notion of |Prendergast| ( |l962[ ) thai 
dust lanes are shocks in the gas hows. She also showed 
that dust lanes are more straight when the bar potential 
is stronger. Due to lack of numerical resolution and/or 
large numerical viscosity, however, early numerical works 
published before the middle of 1990s were unable to cap- 
ture the formation of nuclear rings and spirals accurately. 



Piner et al. ( 1995 ) was the first who focused on the ring 
formation using a grid-based code with negligible numer- 
ical viscosity. Unfortunately, however, their numerical 
results were, at least quantitatively, contaminated by a 
trivial sign e rror in the e valuat ion of the bar fo rces, as 
identifie d by Kim et al. (|2012|). More recently, Ann & 



Thakurj p005| and [Tha 



kur et al. (2009) used SPH sim 



ulations to study gas dynamics in the central regions of 
barred galaxies, but they were unable to resolve nuclear 
spirals because most particles in their simulations are 
captured into a ring, leaving only a small number of par- 
ticles inside the r ing. 

here after Paper I) 
1995) and revisited 



Very recently, Kim et al. (2012 



corrected the error m Piner et al. 



the issue of the bar-substructure formation using high- 
resolution simulations by varying the gas sound speed c s 
and the mass of a central black hole (BH). Paper I found 
that the corrected bar force naturally allows the develop- 
ment of nuc l ear sp irals inside a ring, which was unseen in 



Piner et al. (1995). In addition, Paper I showed that the 



ring size and the mass inflow rate M toward the galaxy 
center depend rather sensitively on c s . In models with 
small c s , the rings are narrow and located away from the 
center. This prevents further inflows of the gas to the 
central regions, making M smaller than 10 -4 M yr~ x . 
On the the other hand, models with large c s have a small 
and broad ring, resulting in M £ 10~ 2 M Q yr~ 1 . Paper I 



furthe r confirmed that the prediction of |Buta fc Combes 
( 1996 ) that the shape of nuclear spirals depend on the 
sign of the d(Jl — K,/2)/dr curve such that they are trail- 
ing (leading) where d(il — n/2)/dr is negative (positive). 
This suggests that nuclear spirals, if exist, are likely to 
be trailing in galaxies with supermassive BHs. 

Since the models studied in Paper I were restricted to 
a particular set of the bar parameters, namely, with the 
bar mass ~ 30% of the spheroidal component (i.e., bulge 
plus bar) and the aspect ratio 1Z = 2.5, they were unable 
to study the variations in the size and shape of bar sub- 
structures with differing bar strength. In addition, the 
models in Paper I are all non-self-gravitating, so that 
the effects of self-gravity on bar substructures h ave yet 
to be explored. Although Athanassoula (| 1992b| consid- 
ered diverse bar models with different aspect ratio and 
bar quadrupole moment Q m , her models were unable to 
resolve nuclea r rings and spirals due to insuf ficient reso- 
lution. While iRegan & Teubenl (120031 120041) also ran a 



large number of numerical simulations with varying Q rn 
and 1Z, their numerica l results were comp romised by the 
incorrect bar forces of Piner et al. (1995), as mentioned 
above. 

In this paper, we present the results of a series of hy- 
drodynamic simulations to investigate the properties of 
bar substructures that form. This work extends Paper 
I in two ways: (1) by considering various bar models 
with differing mass and aspect ratio and (2) by including 
gaseous self-gravity, while fixing the sound speed and the 
BH mass. We measure the curvatures of dust lanes, sizes 
of nuclear ring, shapes of nuclear spirals, and mass inflow 
rates in simulations, and study their relationships with 
the bar parameters. Our main objective is to explore the 
parametric dependence of the properties of bar substruc- 
tures on the bar strength and self-gravity. This will allow 
us to provide physical explanations for the formation of 
bar substructures especially for nuclear rings and spi- 
rals, which was previously uncertain. We also compare 
our numerical results with observations of barred-spiral 
galaxies that are currently available. 

Our models with varying bar strength are useful to ad- 
dress an important question as to what determines the 
position of nuclear rings. Observations indicates that nu- 
clear rings are typically located near the inner Lindblad 
resonance (ILR) when there is a single ILR, or between 
the ILRs when there are two ILRs (e.g., Combes fc Gerin 
1985||Knapen et al.|1995"l|Comer6n et al |2010| ). This has 
often been interpreted as an indication that nuclear rings 
form as a conseque nce of the resonant interactions of the 
gas with the ILRs ([ Combes 1996] |Buta fc Combes||1996 
see also Regan fc 1'eubcn 2003 f. This idea was moti- 
vated by the notion that the bar torque changes its sign 
whenever crossing each ILR such that in the case of a 
single ILR, for example, the gas inside (outside) the ILR 
receives a positive (negative) torque and thus moves out- 
ward (inward). This idea of the resonance-driven ring 
formation requires that the bar torque dominates ther- 
mal and ram pressures of the gas in governing the gas 
motions both inside and outside the ILR. It is clear that 
the bar torque dominates outside the ILR, inducing dust- 
lane shocks and initiating radial gas inflows. However, it 
is uncertain if that is also the case inside the ILR. The bar 
potential becomes increasingly axisymmetric toward the 
galaxy center, while the inflowing gas usually has large 
momentum and is thus unlikely to stall at the ILR. By 
measuring the ring positions and comparing them with 
the ILR radii in our models, we directly test whether 
the ILR is really responsible for the formation of nuclear 
rings. 

This paper is organized as follows: In Section 2, we 
describe our model parameters and numerical methods. 
In Section 3, we define and evaluate the bar strength Qh 
of our galaxy models as well as the radius r(Q(,) where 
the maximum bar torque occurs. We provide algebraic 
expressions for Q\, and r(Qb) for various bar models. In 
Section 4, we quantify the properties of bar substructures 
and explore their correlations with the bar strength. We 
also measure the mass inflows rates. Finally, in Section 5, 
we summarize our results and discuss them in comparison 
with observations. 



2. MODEL 



Gaseous Structures in Barred Galaxies 
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We study gas responses to an imposed non- 
axisymmetric bar potential using hydrodynamic simu- 
lations with and without self-gravity, focusing on the ef- 
fects of the bar strength on gaseous structures that form. 
We consider an initially-uniform, rotating gaseous disk 
with surface density S — 10 M Q pc~ 2 . We assume that 
the disk is infinitesimally thin, isothermal, and unmag- 
netized. The simulations setups and numerical methods 
are identical to those in Paper I, except that some models 
in the present work include self-gravitational potential of 
the gas. Here we briefly summarize our numerical models 
and highlight the differences between the current models 
and those studied in Paper I. 

The gaseous disk is placed under the external gravi- 
tational potential $ oxt consisting of four components: a 
stellar disk, a stellar bulge, a non-axisymmetric stellar 
bar, and a central BH. The stellar disk is modeled by 
a Kuzmin-Toomrc profile, the bulge by a modified Hub- 
ble profile, and the BH by a Plummer pr ofile with mass 
AIbh- For the bar potential $bar, we use 
prolate spheroids with volume density 



Ferrers (1887) 



P 



_ J Pba 



1- 

0, 



9 



for g < 1, 
elsewhere. 



(1) 



where pbar is the central density, g 2 = y 2 /a 2 + (x 2 + 
z 2 )/b 2 , and a and b (< a) denote the semimajor and 
scmiminor axes of the bar, respectively. The index n 
represents the degree of central density concentration. 
Appendix presents an analytic expression of $bar at z = 
and g < 1 for n = 1. 

The total mass of the bar is given by Mh ar = 
2 2n+3 TTab 2 p haT T( n + l)T(n + 2) /T(2n + 4) (e.g., |Athana~ 
soula et al.|2009 ) , which we control using a dimensionless 
parameter 



./"bar 



Mi 



bar 



M bar + M bu l ' 



(2) 



where Mb u i denotes the bulge mass inside r = 10 kpc. In 
our models, the total mass of the spheroidal component 
within 10 kpc is fixed to Af bar + M bu i = 4.87 x 10 10 M Q , 
so that larger /b ar implies a less massive bulge. Another 
independent parameter that characterizes the bar poten- 
tial is the aspect ratio TZ = a/b. One may also use the 
bar quadrupole moment Q m — M b ar (a 2 — b 2 )/(5 + 2n) 



as a measure of the bar mass (e.g., Athanassoula] 1992b 
Regan &: Teuben||2003| |2004| ), although we prefer / bar 
since Q m depends also on TZ. The rotation curve in the 
bar regions is slightly changed by varying /bar and TZ, 
but its fiat part has v = 200 km s _1 regardless of /b ar 
and TZ, corresponding to normal disk galaxies. 

Since the parameter space is very large, Paper I fixed 
the bar parameters to /bar = 0.3, a = 5 kpc, and 
6 = 2 kpc, and varied c s from 5 to 20 km s" 1 and Mbh 
from to 4 x 10 s M to explore how gas flows change 
with the gas sound speed and the BH mass. In this work, 
we instead fix c s = 10 km s _1 and Mbh = 4 x 10 7 M , 
and consider 40 models with differing /b ar = 0.08 to 0.6 
and TZ — 1.5 to 3.5. Since the bar semimajor axis is fixed 
to a = 5 kpc in all of our models, the latter is equivalent 
to varying b from 3.33 to 1.43 kpc. Other parameters 
such as the bar pattern speed fib = 33 km s _1 kpc -1 
and the total disk mass A/disk = 2.2 x 10 11 M Q remain 
the same as in Paper I. Columns (1) - (3) of Tables 



[T] and [2] list the model names and parameters for non- 
self-gravitating and self-gravitating models, respectively. 
The postfixes N and G stand for no self-gravity and with 
self-gravity, respectively. In what follows, a model name 
without any postfix refers to both self-gravitating and 
non-self-gravitating models. Note that Model M30R25N 
is identical to Model csl0bh7 studied in Paper I except 
for the size of the inner boundary (see below). 

With fit, = 33 km s _1 kpc" , all the models have the 
corotation resonance at tcr = 6 kpc, independent of /bar 
and TZ. The presence of a central BH makes the rotation 
curve rise steeply toward the center as oc (MbhA) 1 / 2 , 
giving rise to a single ILR where fib = fi — k/2. Here, 
ft and k refer to the angular and epicycle frequencies, 
respectively, averaged between on the major and minor 
axes of the bar. Unlike rcR, the ILR position hlr de- 
pends slightly on /b ar and TZ. Column (4) of Tables [l] 
and [2] gives hlr m each model. Since the bulge poten- 
tial is more centrally concentrated than the bar potential, 
models with larger /b ar tend to have a smaller rotational 
velocity in the central regions and hence smaller ?'ilr,. 
For fixed /bar, a larger value of TZ corresponds to larger 
Pbar and thus slightly larger ri LR . 

For models in wh ich self-gravity is included , we employ 
the method used by |Shetty fc Ostriker ( 2008 1 to calculate 
the gravitational potential trom the perturbed gas sur- 
face d ensity, S — SqQ This method is based on Kalnajs 



(1971 1 's scheme that is efficient on a logarithmically 



spaced radial grid. As a softening parameter, we take 
H/r = 0.1 that allows for non-zero thickness of the disk 
in the potential calculation. In our models, the gaseous 
disk initially has a Toomre stability parameter of 



2.1 



vrGE 



10 kpc 



10 M pc- 



CS) 

so that it is gravitationally stable. However, nuclear rings 
that form at around r ~ 1 kpc would achieve sufficient 
density to be unstable when the bar is strong. 

We solve the basic equations of ideal hydrodynamics 
using the CMHOG code in a frame corotating with the 
bar (Paper I). CMHOG is third-order accurate in spac e 
and has very little numerical diffusion ( Piner et al.|1995 1. 
We resolve the central regions with high accuracy by 
setting up a logarithmically-spaced cylindrical grid ex- 
tending from 0.04 kpc at the inner boundary to 16 kpc 
at the outer boundary. The size of the inner bound- 
ary taken in the present models is twice larger than that 
in Paper I: we checked that this makes negligible differ- 
ences in the properties of bar substructures except for the 
mass inflow rate across the inner boundary that becomes 
larger, by ~ 10% on average, in the current models. The 
number of zones in our models is 1024 in the radial di- 
rection and 535 in the azimuthal direction that covers 
the half-plane. The corresponding spatial resolution is 
Ar = 0.23, 5.86, and 93.8 pc at the inner boundary, 
r = 1 kpc, and the outer boundary, respectively. We 
adopt the outflow and continuous boundary conditions 
at the inner and outer boundaries, respectively, while 

1 In our self-gravitating models, gravity of the initial gaseous 
disk with surface density So is neglected in order to make the 
initial rotation curve identical to that in the non-self-gravitating 
counterparts. 
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Table 1 

Model Parameters and Curvatures of Dust Lanes in Non-self-gravitating Models 



Model 


/bar 


11 


r ILR 


Qt 


r(Qb) 


Oil 




Q 2 




Aa 








(kpc) 




(kpc) 


(< 


eg 


) 


(deg) 


(deg) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 




\°) 




(9) 


M08R15N 


0.08 


1.5 


2.2 


0.02 


3.60 


70.5 


± 


0.5 


35.8 ± 


0.9 


64.2 ± 2.5 


M08R20N 


0.08 


2.0 


2.2 


0.05 


3.30 


74.4 


± 


5.8 


35.8 ± 


2.9 


63.4 ± 12.6 


M08R25N 


0.08 


2.5 


2.2 


0.07 


3.12 


74.1 


± 


2.6 


44.2 ± 


2.3 


46.3 ± 3.4 


M08R30N 


0.08 


3.0 


2.2 


0.09 


3.04 


72.8 


± 


3.0 


53.3 ± 


4.0 


30.1 ± 6.6 


M08R35N 


0.08 


3.5 


2.2 


0.12 


2.98 


74.3 


± 


1.7 


59.6 ± 


5.2 


22.2 ± 6.7 


M15R15N 


0.15 


1.5 


2.1 


0.04 


3.58 


75.1 


± 


1.4 


34.5 ± 


1.6 


81.6 ± 2.9 


M15R20N 


0.15 


2.0 


2.1 


0.08 


3.28 


71.7 


± 


3.2 


49.6 ± 


4.3 


36.3 ± 9.0 


M15R25N 


0.15 


2.5 


2.2 


0.13 


3.12 


75.3 


± 


2.5 


63.4 ± 


2.3 


19.6 ± 4.2 


M15R30N 


0.15 


3.0 


2.2 


0.17 


3.04 


73.4 


± 


2.5 


68.3 ± 


2.4 


8.2 ± 4.4 


M15R35N 


0.15 


3.5 


2.2 


0.21 


2.98 


74.3 


± 


5.1 


71.1 ± 


5.4 


5.0 ± 3.6 


M30R15N 


0.30 


1.5 


1.8 


0.08 


3.56 


70.7 


± 


3.4 


38.9 ± 


3.9 


63.4 ± 10.0 


M30R20N 


0.30 


2.0 


1.8 


0.16 


3.26 


78.8 


± 


3.1 


60.6 ± 


4.5 


33.4 ± 9.6 


M30R25N 


0.30 


2.5 


2.1 


0.23 


3.10 


79.8 


± 


1.5 


68.5 ± 


1.4 


16.8 ± 2.8 


M30R30N 


0.30 


3.0 


2.1 


0.31 


3.02 


78.4 


± 


2.8 


75.1 ± 


0.9 


4.7 ± 3.9 


M30R35N 


0.30 


3.5 


2.1 


0.39 


2.90 


80.8 


± 


0.3 


79.3 ± 


0.5 


2.2 ± 0.2 


M60R15N 


0.60 


1.5 


1.4 


0.14 


3.42 


81.0 


± 


3.3 


52.6 ± 


11.1 


53.4 ± 23.1 


M60R20N 


0.60 


2.0 


1.3 


0.28 


3.18 


79.2 


± 


3.1 


70.7 ± 


1.1 


12.8 ± 5.7 


M60R25N 


0.60 


2.5 


1.3 


0.41 


3.02 


79.3 


± 


3.4 


74.8 ± 


1.9 


6.9 ± 5.0 


M60R30N 


0.60 


3.0 


1.9 


0.54 


2.92 














M60R35N 


0.60 


3.5 


1.9 


0.67 


2.88 















Note. — /tax is the fraction of the bar mass relative to the spheroidal component (bar 
plus bulge); TZ is the ratio of the bar semimajor axis to the semiminor axis; t-tlr is the 
radius of the ILR; Q;, is the bar strength defined by equation El; r(Qb) is the radius of 
the maximum bar torque; c?i and 0.2 are the tangent angles to a aust-lane segment at the 
inner and outer ends relative to the x-axis, respectively; Ac? is the dust lane curvature 
defined by equation J7j). 



Table 2 

Model Parameters and Curvatures of Dust Lanes in Self-gravitating Models 



Model 


/bar 


K 


Hlr 


Qb 


r(Q b ) 


Oil 




«2 




Aa 








(kpc) 




(kpc) 


(deg) 


(deg 


) 


(deg) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 




(8) 




(9) 


M08R15G 


0.08 


1.5 


2.2 


0.02 


3.60 


69.7 ± 


1.5 


37.1 ± 


1.5 


60.2 ± 1.8 


M08R20G 


0.08 


2.0 


2.2 


0.05 


3.30 


72.1 ± 


7.6 


35.8 ± 


2.6 


57.1 ± 15.0 


M08R25G 


0.08 


2.5 


2.2 


0.07 


3.12 


70.7 ± 


3.8 


47.2 ± 


10.1 


39.9 ± 20.8 


M08R30G 


0.08 


3.0 


2.2 


0.09 


3.04 


69.5 ± 


2.1 


52.6 ± 


2.7 


25.5 ± 5.4 


M08R35G 


0.08 


3.5 


2.2 


0.12 


2.98 


73.6 ± 


4.5 


61.4 ± 


6.7 


28.9 ± 5.3 


M15R15G 


0.15 


1.5 


2.1 


0.04 


3.58 


74.2 ± 


4.0 


39.1 ± 


1.5 


74.7 ± 13.9 


M15R20G 


0.15 


2.0 


2.1 


0.08 


3.28 


71.0 ± 


1.1 


50.3 ± 


1.7 


32.8 ± 13.2 


M15R25G 


0.15 


2.5 


2.2 


0.13 


3.12 


75.0 ± 


3.5 


65.6 ± 


2.3 


16.1 ± 4.4 


M15R30G 


0.15 


3.0 


2.2 


0.17 


3.04 


72.9 ± 


0.9 


69.0 ± 


3.0 


6.3 ± 4.1 


M15R35G 


0.15 


3.5 


2.2 


0.21 


2.98 


74.7 ± 


1.6 


69.8 ± 


3.5 


8.0 ± 6.8 


M30R15G 


0.30 


1.5 


1.8 


0.08 


3.56 


74.5 ± 


4.2 


42.7 ± 


3.6 


65.7 ± 19.8 


M30R20G 


0.30 


2.0 


1.8 


0.16 


3.26 


75.1 ± 


5.9 


60.6 ± 


1.9 


25.2 ± 11.1 


M30R25G 


0.30 


2.5 


2.1 


0.23 


3.10 


78.5 ± 


2.6 


69.7 ± 


2.8 


18.2 ± 10.6 


M30R30G 


0.30 


3.0 


2.1 


0.31 


3.02 


78.2 ± 


3.2 


75.2 ± 


5.0 


6.4 ± 6.0 


M30R35G 


0.30 


3.5 


2.1 


0.39 


2.90 


79.6 ± 


0.1 


79.7 ± 


0.5 


2.7 ± 1.2 


M60R15G 


0.60 


1.5 


1.4 


0.14 


3.42 


80.2 ± 


3.1 


57.1 ± 


5.7 


64.2 ± 24.0 


M60R20G 


0.60 


2.0 


1.3 


0.28 


3.18 


79.1 ± 


1.1 


69.3 ± 


8.9 


22.8 ± 16.2 


M60R25G 


0.60 


2.5 


1.3 


0.41 


3.02 


76.6 ± 


3.2 


72.7 ± 


1.5 


9.2 ± 9.8 


M60R30G 


0.60 


3.0 


1.9 


0.54 


2.92 












M60R35G 


0.60 


3.5 


1.9 


0.67 


2.88 













Note. — See Table [T] for the definitions of the various symbols 



taking the periodic conditions at the azimuthal bound- 
aries. The bar potential is slowly turned on over one bar 
revolution time 2tt/ fib = 186 Myr to minimize transients 
in the flows caused by its sudden introduction. All the 
models are run until 1 Gyr. 

3. BAR STRENGTH 

One of the key parameters that govern the gas dynam- 
ics in barred galaxies is the bar strength which depends 
on both /bar and 1Z. In this section, we evaluate the 



bar strength of our numerical models that will be used 
to analyze the properties of bar substructures in Section 
4. We also provide fitting formulae for the bar strength 
and the position of the maximum bar torque for future 
purposes. 

It has often been customary to measure the bar 
strength using the dimensionless parameter Qb defined 
by 

^ _ F T {r,4 



Fr(t) 



(4) 
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[T] and [2] list Qb and r(Qb) of our numerical models. 
By varying /bar, TZ, and n, we find that 



Figure 1. Dependence of the bar strength Q;, and the maximum- 
torque radius r(Q(,) on the bar aspect ratio TZ = a/b for /i, ar = 0.3. 
The solid curves plot the numerical results for our galaxy models 
with n = 0,1,2 Ferrers bars, while the dashed lines are for the 
analytic results assuming flat rotation. 

where Fx = — {d^^/dcfy/r is the tangential force 
due to the non-axisymmetric bar potential and Fr = 

— (l/27r) Jq* d$ ext /drd(j) is the azimuthally-averaged ra- 
dial fOTce (ej^_JCombeF&^^ 

et al.||2UU41 ILaurikainen et al.||2004 

aTl |2009 



Salo||2002| 
20061 JPeep 



Block 



Marti 



et 



| |2006| IGomCTon - 
2010p Physically, Qb corresponds to the maximum bar 



torque applied to a gaseous materia l in orbital motion 
relative to its specific kinetic energy. Block et al.| ( 2004 ) 
found that Qb is deeply related to the class of barred- 
spiral galaxies such that SA galaxies have Qb £ 0.1, while 
Qb <L 0.15 for SB galaxies. 

Since the bar torque is stronger for a more massive 
and elongated bar, Qb should be an increasing function 
of /bar an d TZ. It is straightforward to calculate Qb values 
of our galaxy models as well as the radius r(Qb) where 
\Ft\/Fr is maximized. Figure [l] plots as solid lines the 
resulting Qb and r(Q&) as functions of TZ for the Ferrers 
bar with /b ar = 0.3 and n = 0,1,2. Note that when 
/bar = 0.3, Qb for n = 2 is not much different from the 
case with n = 1, while r(Qb) becomes smaller by ~ 16% 
as n increases from 1 to 2. Also plotted in Figure [T] 
as dashed lines are the analytic expressions derived in 
Appendix for the Ferrers bars with n — 1 assuming flat 
rotation, which agrees with the numerical results within 
less than 15%. The difference between the numerical and 
analytic results for n — 1 is of course due to the fact that 
the rotation curves in our models are not strictly flat but 
decrease slowly with r in the regions of the maximum bar 
torque, tending to increase r(Qb) compared to the case 
with exactly flat rotation. Columns (5) and (6) of Tables 

2 The non-axisymmetric torque is sometimes rep resented by Q„ 
that includes a contributio n from spiral arms (e.g., [Durbala et al.| 
|2009|[Comer6n et al.|2010| . In our models, Qf, = Q g since no spiral 
arm is considered. 



•58/£f (a/6-1), forn = 0, 
.44/&f(a/&-l), forn = l, 



0., 

0.38/^ (a/6-1), forn = 2, 



(5) 



and 



1.024 
0.934 
0.817 



0.161(a/6) + 0.024(a/6) 2 , for n = 0, 
Q.196(o/6) + 0.028(a/6) 2 , for n = 1, 
0.197(a/6) + 0.031(a/6) 2 , for n = 2, 

(6) 

are good fits, within 5%, to the numerical results. While 
Qb is linearly proportional to 1Z, r{Qb) depends very 
weakly on 1Z and is independent of /bar- The reason 
for Qb varying less steeply than /^ a ° is that large /b ar 
increases the rotational velocity slightly near the regions 
where the bar torque achieves its maximum. 

4. SIMULATION RESULTS 

Paper I presented gas dynamical evolution of Model 
M30R25N with Q b = 0.23' (/bar = 0.3 and TZ = 2.5) in 
detail. Evolution of other models with differing Qb are 
qualitatively similar except that dust lanes and/or nu- 
clear rings are absent if /b ar and/or 1Z are too large. In 
this section, we first briefly describe the main evolution- 
ary features and the conditions for the existence of dust 
lanes and nuclear rings, and then present detailed anal- 
yses of the properties of bar substructures formed in our 
models. 

4.1. Overall Morphology 

As the non-axisymmetric bar potential is slowly turned 
on, initially-circular gas orbits are perturbed due to the 
bar torque, creating overdense ridges at the downstream 
side of the bar major axis. As the amplitude of the 
bar potential grows further, the ridges soon develop into 
dust-lane shocks. Dust-lane shock s are found unstable to 
the wiggl e instability identified b y Wada & Koda ( 2004[ ) 
(see also Kim fc Ostriker||2006[ ) and form clumps with 
high vorticty (see Paper I for more detailed description). 
Gas loses angular momentum at the shocks, flows radially 
inward, and forms a nuclear ring at the position where 
the centrifugal force balances the external gravity. Sonic 
perturbations launched from the ring propagate inward 
and excite m = 2 nuclear spirals in the nuclear regions. 
With Mbh = 7x 10 7 M , nuclear spirals that persist in 
our models are all trailing. 

Figures [2] and [3] display snapshots of gas surface den- 
sity in logarithmic scale at t = 0.3 Gyr from all non- 
self-gravitating and self-gravitating models, respectively. 
Figure [4] plots density distributions at t — 0.8 Gyr 
for both non-self-gravitating (two upper rows) and self- 
gravitating (two bottom rows) models with 7Z = 2.0 
and 3.0. The inner ±5 kpc regions are shown. The 
bar is oriented vertically along the y-axis. The solid 
lines in each panel plot xi-orbits that cut the s-axis at 
x c = 0.8,1.4,2.0, and 3.6 kpc under the total external 
potential. It is apparent that gas responses inside the 
outermost X!-orbit with x c = 3.6 kpc are quite dramatic, 
while the outer reg ions are relatively unperturbed (e.g., 
Kim fc StonepOlip . 

At t = 3 Gyr, all the models with Q b < 0.5 
(/bar < 0.3 or 72. < 2.5) contain well-defined dust lanes 



G 



Kim et al. 




Figure 2. Logarithm of the gas surface density at t = 0.3 Gyr for non-self-gravitating models. Each panel shows the inner ±5 kpc regions 
where the bar is oriented vertically along the y-axis. Each row corresponds to models with 1Z = 1.5, 2.0, 2.5, 3.0, and 3.5 from top to 
bottom, while each column is for models with /t, ar = 0.08, 0.15, 0.3, and 0.6 from left to right. Solid lines in each panel draw xi-orbits 
that cut the rr-axis at x c = 0.8, 1.4, 2.0, and 3.6 kpc. Color bar labels log(S/So). 
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Figure 4. Logarithm of the gas surface density at t = 0.8 Gyr for all models with TZ = 2.0 and 3.0. The two upper rows plot the 
non-self-gravitating models, while the two bottom rows are for the self-gravitating models. Solid lines in each panel draw xi-orbits that 
cut the x-axis at x c = 0.8, 1.4, 2.0, and 3.6 kpc. Color bar labels log(S/Eo)- 



that follow an xi-orbit reasonably well, although they 
become weaker aXt — 0.8 Gyr. On the other hand, Mod- 
els with Qb > 0.5 do not possess features that resemble 
observed dust lanes. In these high-Qb models, the bar 
torque is so strong that the gas inside the outermost X\- 
orbit loses most of its initial angular momentum even 
in the developing stage of dust lanes (t ~ 0.13 Gyr for 
Model M60R35N). With the shocked gas rapidly lost to 
the nuclear regions and subsequently to the inner bound- 
ary, there remains no material available to support the 
dust lanes in the bar regions. In these models, the bar 



regions are occupied by gas blobs in a filamentary shape 
that are repeatedly stretched and folded on the course of 
orbital motions along xi-orbits about the center. Com- 
parison between Figures [2] and [3] shows that the shape of 
dust lanes is not much affectedby self-gravity. 

Figures [2] also shows that conventional a^-type rings, 
that is, rings whose shape is similar to a^-orbits, exist 
when /b ar < 0.3 or 1Z < 1.5, in non-self-gravitating mod- 
els. Rings are quite clumpy even in non-self-gravitating 
models due to addition of clumps produced by the wig- 
gle instability at the dust-lane shocks. Nuclear rings in 
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Figure 5. Temporal evolution of (a) the peak gas surface density 
Speak an d (b) the ^-coordinate Xgcak of the dust lanes at r = 
1.5 kpc for Models M08R15, M15R20, and M30R25. The solid 
and dotted lines are for the self-gravitating and non-self-gravitating 
models, respectively. For Models M30R25N and M30R25G, only 
the results up to t = 0.54 Gyr are shown, after which dust lanes 
are too weak to be defined well. 



Models M60R20N and M60R25N are highly inclined rela- 
tive to the x-axis at t = 0.3 Gyr and gradually precess to 
align their long axes to the bar major axis, eventually be- 
coming xi-type rings, although they become very weak 
at late time (e.g., Fig. |4[i). In Models M60R30N and 
M60R35N, orbits of gas m the nuclear regions are highly 
transient to be considered as rings. In self-gravitating 
models, on the other hand, rings that form are all of the 
x 2 type (Fig. [3}. In Models M60R30G and M60R35G, 
rings at t = O.oGyr are so small that they become weaker 
with time by losing mass directly to the central hole, dis- 
appearing at t ~ 0.75 Gyr. For models that form an 
22-type ring, rings are larger in size in self-gravitating 
models than in non-self-gravitating models (Fig. |4| . We 
will explain how gas orbits are affec ted by the bar poten- 
tial and self-gravity in Section 4.3 below. 

As Qb decreases, angular momentum loss at the shocks 
becomes smaller, resulting in a larger nuclear ring. Dust 
lanes are correspondingly located further downstream 
from the bar major axis since their inner ends are always 
attached to a nuclear ring. While dust lanes become 
weaker with time after the peak strength is attained, nu- 
clear rings remain strong until the end of the runs (e.g., 
Fig. [4 l, suggesting that they are long-lived features (e.g . 



Allard et al. 120061 ISarzi et al.|2007 Comeron etal.|2010). 



1 



In Models M08R15N and M08R15G with Q b = 0.02, the 
bar torque is so weak that the amount of gas moving in 
along the dust lanes is not sufficient to form an appre- 
ciable ring until the end of the runs. 

The presence of nuclear spirals is deeply related to the 
shape and size of a nuclear ring. There is no chance to 
possess nuclear spirals in models with an eccentric x\- 
type ring. Even in models with an a^-type ring, nuclear 
spirals are absent if a nuclear ring is too small (as in Mod- 



els M30R35N, M30R35G, M60R20G, and M60R25G) to 
provide enough space for coherent structures to grow. In 
Models M08R15N and M08R15G, on the other hand, the 
inflowing gas excites trailing spiral waves at r ~ 1 kpc 
that propagate radially inward and turn into nuclear spi- 
rals. We defer the more detailed disc ussio n on temporal 
evolution of nuclear spirals to Section |4~4| 

4.2. Dust Lanes 

To quantify the strength and displacement of dust 
lanes from the bar major axis, we confine to the regions 
with r = 1.5 kpc and measure the ^-coordinate, x pea k, 
of the position where gas surface density is maximized 
at each time. Figure [5] plots the temporal changes of the 
peak surface density £ pca k as well as £ pe ak at r = 1.5 kpc 
in Models M08R15, M15R20, and M30R25. The solid 
and dotted lines are for the self-gravitating and non- 
self-gravitating models, respectively. In Model M30R25, 
Sp ea k rises with time and peaks at t ~ 0.13 — 0.15 Gyr 
after which it decreases with time as the bar regions be- 
come increasingly evacuated due to gas infalls to the ring. 
At early time when the bar potential is weak, dense ridges 
form at x pC ak ~ —2.5 kpc which move closer to the bar 
major axis as they turn to dust-lane shocks and stay at 
£peak ~ — 1 kpc. In these models with Qb = 0.25, the 
decay of the dust lanes is relatively rapid that they be- 
come almost invisible after t = 0.54 Gyr. The temporal 
behaviors of £ pea k and £ pea k are qualitatively the same 
for other models with different bar strength, although 
the dust lanes tend to decay more slowly and locate far- 
ther from the bar major axis as Qb decreases. Note that 
while the strength of dust lanes varies with time con- 
siderably, their location and shapes do not change much 
after t = 0.25 Gyr. Note also that self-gravity does not 
make much changes in £ pca k and 2 pea k- This is because 
dust lanes with typical density S pea k ~ 100 M Q pc~ 2 at 
r = 1.5 kpc have Qt ~ 1.4, larger than unity (e.g., Eq. 
[3]). In addition, dust lanes are known to possess strong 
velocity shear, about 10 times larger than the velocity 
shear in disks at large (Paper I), which tends to reduce 
t he mass-collecting effec t of self-gravity in dust lanes. 



Comeron et al. (20091 defined the dimensionless cur- 
vature of a dust-lane segment as 



Aa 



a 2 



r(Q b ), 



(7) 



where a\ }2 and r 12 indicate the tangent angle to, and 
the position vector of, the inner and outer ends of the 
se gment, respectively. T his extends the definition of Aa 



of the host bar. 



Knapen et al. ( 2002 ) to take i nto account the length 



(2009) 



lb measure Aa, Comeron et al 
selected a constant-curvature part of a dust lane in each 
galaxy, with the determination of both ends relying on 
visual inspection. In order to obtain Aa unambigu- 
ously from our numerical models, we instead confine to 
a dust-lane segment bounded by the fixed position an- 
gles 4>i t2 = 135°, 100° (or, 4> h2 = -45°, -80°) from the 
positive x-axis, and measure ct\ 2 and ri 2 at both ends. 
Columns (7)-(9) of Tables [11 and B list the mean values 
(together with standard deviations) of a,\, a 2 , and Aa 
averaged over t = 0.25-0.35 Gyr for models that possess 
dust lanes. While a 2 varies rather sensitively with Qb, 
ai ~ 70°-80° does not change much, indicating that Aa 
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Figure 6. Curvatures Act of dust lanes as functions of (a) Qt, alone and (b) a linear combination of Qf, and Ti = a/b for all models that 
posses dust lanes. Various symbols give the mean values of Aa averaged over t = 0.25-0.35 Gyr, while the errorbars indicate the standard 
deviations. In both panels, filled and open symbols are for self-gravitating and non-self-gravitating models, respectively. The dotted line 
in (b) draws the best fit, expressed in equation JsJ. 
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Figure 7. Dependence of xi-orbits with x c = 0.6, 1.0, 1.4, 2.0 kpc (a) on Q;, when 1Z = 2.0 is fixed, and (b) on 1Z when 
Note that the parts a;i-orbits in a quadrant are more curved as x c increases or TZ decreases. 

is determined primarily by the tangent angle at the outer 
end of dust lanes. 

FigurelB^, plots Aa measured from our numerical mod- 
els as a function of Qt,. Filled and open symbols are for 
self-gravitating and non-self-gravitating models, respec- 
tively. These agree within error bars, indicating that self- 
gravity does not affect the shape of dust lanes much. In 
general, Aa decr eases with increasing Q b, consist ent with 
the prediction of Athanassoula ( 1992b ) (see also Knapen 



0.2 is fixed. 



et al. 2002). Note however that there is a considerable 
scatter in Aa for given Qt- Comeron et al. (2009) found 
that the spread in the data points can be much reduced 
by fitting them using a linear combination of Aa and 1Z. 
We follow the same procedure and find the best fit 



Q b + 0.10ft = 0.87 - 0.37 log Aa, 



(8) 



for all models, which is plotted as a dotted line in Figure 
[(jjx The chi-square measure of equation (|8| is about 3 
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Figure 8. Snapshots of logarithm of surface density at t = 0.3 Gyr in the inner 2 kpc regions of non-self-gravitating models (a) M08R20N, 
(b) M30R20N, and (c) M60R20N, and self-gravitating models (d) M08R20G, (e) M30R20G, and (f) M60R20G. The solid lines draw x 2 - 
orbits that cut the x-axis at x c = 0.8,1.2,1.6 kpc in (a), (b), (d), and (e), xi-orbits with x c = 0.2,0.4,0.6 kpc in (c), and X2-orbits with 
x c = 0.8, 1.2 kpc in (f). Color bar labels log(E/E ). 



times smaller than that of the best fit using only Qb- 
Equation ^ implies that more elongated bars lead to 
more straight dust lanes. 
The behavior of Aa upon Qb and TZ given in equation 
can be qualitatively understood from the facts that 
lust lanes roughly follow xi-orbits and that they move 
closer to the bar major axis as Qb increases. Figure [7] 
illustrates the changes in the shapes of xi-orbits with 
x c = 0.6,1.0,1.4,2.0 kpc as Qb or TZ varies. For fixed 
1Z, the parts of xi-orbits, say, in the second or fourth 
quadrant where dust lanes are located, are clearly more 
curved as x c increases, explaining larger Aa in models 
with smaller Qb- For fixed Qb, si-orbits are more acute 
near the y-axis under a more elongated bar (i.e., larger 
TZ), while they are largely independent of TZ near the x- 
axis. This makes ot\ almost unchanged with TZ, while 
causing a2 to increase with TZ, resulting in smaller Aa 
for larger TZ. 

4.3. Nuclear Ring 

A conventional wisdom about nuclear rings is that they 
form near the ILRs as a result of resonant interactions of 
gas with the backgrou nd potential (e.g., Combes 1996 
Buta & Combes 1996). In this subsection, we use our 



numerical models to study what controls the ring for- 



mation and its size when there is a single ILR. Paper 
I showed that the mass of a central BH (and thus the 
number and locations of ILRs) does not much affect the 
physical properties of nuclear rings that form, so that 
the results presented below hold also for models with two 
ILRs. We first describe the results of non-self-gravitating 
models and then discuss the effect of self-gravity on nu- 
clear rings. 

4.3.1. Non-self-gravttating Models 

Inspection of Figure [2] reveals that the ring shape and 
its position depend on /bar and/or TZ. To show this more 
clearly, the upper panels of Figure [8] zoom in on the cen- 
tral ±2 kpc regions and plots density distributions in 
logarithmic scale at t — 0.3 Gyr together with a few 
x\- or X2-orbits for selected models with TZ = 2.0. It is 
apparent that Models M08R20N and M30R20N have an 
^2-type nuclear ring, while the ring in Model M60R20N 
is inclined in the galaxy plane with respect to the x-axis. 
Rings become smaller with increasing /bar, while Model 
M60R20N has a very eccentric ring. 

To illustrate how nuclear rings form in our models, 
we go back to early time when they were beginning to 
shape. Figure [9] plots the instantaneous streamlines of 
the gas that starts from Point A marked at (x,y) — 
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Figure 9. (a) Instantaneous streamlines of the gas that starts 
from Point A (x,y) = (0,2.5 kpc) in Model M08R20N at t = 
0.25 Gyr (red), Model M30R20N at t = 0.15 Gyr (blue), and 
Model M60R20N at t = 0.12 Gyr (black) when a nuclear ring is 
beginning to form. The dotted circles mark the locations of the 
ILR at n LR = 2.2,1.8,1.3 kpc for Models M08R20N, M30R20N, 
M60R20N, respectively, (b and c) The variations of the azimuthal 
and radial velocities of the gas along the paths shown in (a). The 
ILRs are indicated as vertical lines. In (b), the dashed lines draw 
the equilibrium circular velocity in each model. 

(0, 2.5 kpc) as well as the variations of the rotational 
and azimuthal velocities in the inertial frame along the 
streamlines for Model M08R20N at t = 0.25 Gyr (red), 
Model M30R20N at t = 0.15 Gyr (blue), and Model 
M60R20N at t = 0.12 Gyr (black). In each panel, the 
dotted circles or vertical lines mark the ILR located at 
n LR = 2.2, 1.8, 1.3 kpc for Models M08R20N, M30R20N, 
M60R20N, respectively. Dashed lines in Figure draw 
the equilibrium rotation curve, averaged between on the 
major and minor axes of the bar, in each model. A sim- 
ilar plot is given in Figur e 5 of Paper I for mo dels with 
no BH and in Figure 5 of Kim & Stone (2012) for mag- 
netized models. 
Outside the ILR, the bar provides a negative torque 



for the gas, producing dust-lane shocks. In Model 
M08R20N, the gas that passes through Point A hits the 
shocks at Point B. By losing angular momentum there, 
it starts to move radially in and crosses the ILR at Point 
C. The radial velocity of the inflowing gas at Point C is 
rather large at ~ 63 km s _1 . The associated ram pres- 
sure is larger, by about a factor of 40 and 4, respectively, 
than the thermal pressure and the bar torque at the ILR 
that may try to stop the inflowing motion of the gas. 
Thus, the inflowing gas is not halted at the ILR and con- 
tinues to move inward. At the same time, the gas rotates 
gradually faster at the expense of the gravitational po- 
tential energy. It is at Point D that the gas achieves its 
rotational velocity comparable to the equilibrium value 
and subsequently makes a closed-loop orbit, finally form- 
ing a nuclear ring. The mean radius of the ring in Model 
M08R20N is r ring ~ 1.4 kpc, well inside the ILR. Models 
with a more massive bar forms a ring closer to the center, 
owing to a larger loss of angular momentum. All rings 
form at the position where the centrifugal force balances 
the external gravity. These imply that the formation of 
nuclear rings is not determined by the ILR, but by the 
centrifugal barrier that the gas driven inward by the bar 
torque cannot overcome. 

The upper panels of Figure [10] show the temporal and 
radial variation of the az imut hally- averaged surface den- 
sity for all non-self-gravitating models with 1Z = 2, with 
the horizontal line indicatin g th e ILR in each model. Var- 
ious open symbols in Figure 11 plot the velocities u y , seen 



in the rotating frame with the bar, of the ring gas at y = 
when the ring is beginning to form in the 1Z = 2 models 
without self-gravity. Also plotted are the y-velocities of 
xi-orbits (y Xl ; upper curves) and a^-orbits {y X2 \ lower 
curves) when they cut the x-axis. For fixed 1Z, a model 
with larger /b ar has a shallower external potential <& cxt 
and thus a lower value of y Xl and y X2 at the same x. 
A small solid circle at the tip of each of the lowest two 
curves marks the outermost X2-orbit beyond which no 
a:2-orbit exists in that model. 

As mentioned before, the ring in Model M08R20N with 
/bar = 0.08 begins to form at t ~ 0.25 Gyr. Already at 
this time, the ring material has a velocity very similar to 
that of an a^-orbit at the ring position (open circles in 
Fig. 11 ) . Since the bar torque is quite weak, the angular 
momentum of the gas that is subsequently added to the 
ring is not much different from that of the gas already 
in th e ri ng, keeping r r j ng almost unchanged with time 
(Fig. [lOfc,). As /bar (or Qb) increases, rings form earlier 
and locate closer to the center. For instance, the ring 
in Model M30R20N with f hax = 0.3 begins to form at 
r r i ng ~ 1.2 kpc when t ~ 0.15 Gyr, with u y ~ y X2 <C y Xl 
(open squares in Fig. 111. With a considerable loss of an- 
gular momentum at the dust-lane shocks, the gas newly 
added to the ring has increasingly lower angular momen- 
tum, causing it to shrink with time. In addition, thermal 
pressure perturbations make the orbits of ring material 
deviate from a^-orbits. Some gases at the inner parts 
of the ring even take on si-orbits. This causes the ring 
not only to be distributed more widely but also to shrink 
with time as it loses further angular momentum due to 
supersonic collisions of the gases on x\- and X2-orbits. 
Yet, the decreasing rate of the ring radius is still quite 
Jdt 0.1 Gyr" 1 in Model M15R20N 
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Figure 10. Temporal changes of the azimuthally-averaged surface density, {£), in logarithmic scale for all models with TZ = 2.0. The 
upper row is for the non-self-gravitating models, while the lower low is for the self-gravitating counterparts. The horizontal line in each 
panel represents the location of the ILR. Color bar labels log((£)/Erj). 

and 0.4 Gyr" 1 in Model M30R20N. In these models, 

rings are still of the x-i type (Fig. pb). 

In Model M60R20N with / bar = 076, on the other hand, 
a large bar torque takes the inflowing gas to r r i ng ~ 
0.9 kpc with velocity u y ~ 190 km s , which is just 
in between y Xl = 272 km s _1 and u^, = 109 km s _1 at 
that location (open triangles in Fig. 11 1. Since u y 3> y X2 , 
the inflowing gas is unable to settle on an a^-orbit: the 
ring gas instead follows a hybrid orbit that is inclined 
by about </> r i n g ~ 17° with respect to the x-axis at 
t ~ 0.12 Gyr. With a subsequent addition of low angular- 
momentum gas from outside, the ring shrinks and be- 
comes more eccentric and inclined with time. The ring 
position angle is increased to r i ng ~ 40° at t ~ 0.2 Gyr 
when the ring becomes very eccentric to touch the in- 
ner boundary and starts to lose a significa nt a mount of 
its mass through the inner boundary (Fig. 101). At the 



same time, the inclined ring precesses slowly m the clock- 
wise direction due to the bar torque that tends to align 



the ring parallel to its major axis (see |Kim fc Stone|2012| 
for the precession of a ring in magnetized models). Fig- 
ure [8]; shows that the position angle of the ring in Model 
M60R20N is </> ri n g ~ 65° at t = 0.3 Gyr, which eventually 
becomes 90° at t ~ 0.38 Gyr (i.e., xi-type ring). 

The dominance (or absence) of the x± family of the 
gaseous orbits in our non-self-gravitating simulations ap- 
pears to be determined not only by the amount of an- 
gular momentum loss at dust-lane shocks but also by 
the kinetic energies of x\- and a^-orbits allowed under 
a given external potential. The former is measured by 
Qb, while the latter is affected by 7Z. For example, Fig- 
ure^ shows that Model M30R35N with Q b = 0.39 still 
possesses a well-defined x 2 -type ring. In this model, 
Qb is large enough to decrease the velocity of the ring 
material, while 1Z (together with /bar) sets x\- and X2- 
orbits such that the rotational velocity of the ring mate- 
rial is larger than y X2 only slightly, but much less than 
y Xl , leading to an X2-type ring. With an even larger 
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Figure 11. Distributions of {/-velocities in the rotating frame of 
xi-orbits (upper curves) and X2-orbits (lower curves) at y = as 
functions of x. Various symbols plot the numerical data of the 
nuclear ring at y = when it begins to form at t = 0.25 Gyr for 
Model M08R20 (circles), t = 0.20 Gyr for Model M15R20 (trian- 
gles), t = 0.15 Gyr for Model M30R20, and t = 0.12 Gyr for Model 
M60R20. Open and filled symbols give the results for the non-self- 
gravitating and self-gravitating models, respectively. A small dot 
at the tip of each of the lowest two curves denotes the outermost 
22-orbit in that model. 

value of Q b = 0.67, the ring in Model M60R35N has 
u y ~ 120 km s _1 at r ~ 0.8 kpc when it begins to 
shape at t — 0.08 Gyr. However, this model has the 
outermost a^-orbit at x c — 0.43 kpc, so that the ring 
material has no choice but to follow a^i-orbits. The ring 
subsequently becomes smaller and more eccentric as low 
angular-momentum gas is continuously added. As the 
ring material directly plunges through the inner bound- 
ary, it rapidly decays at a rate d\nr T i n g/dt ~ —15 Gyr -1 . 
In this model, most of the ring gas is accreted to the 
center by t ~ 0.25 Gyr, with a small quantity of gas in 
filamentary shapes moving about the center along X\- 
rather than a^-orbits. 

4.3.2. Effects of Self- gravity 

Nuclear rings form typically at r ~ 0.5 — 1.5 kpc. 
Equation (|3| suggests that a ring could possibly be 
gravitationally unstable if its peak density exceeds ~ 
500 Mq pc _2 p] This happens in most of our self- 
gravitating models except for Models M08R15G and 
M15R15G where the bar torque with Qb ,$ 0.04 is too 
weak to cause significant gas inflows to the central re- 
gions. If a ring were uniform in models with Qb £ 0.05, 
it would have been prone to gravitational fragmenta- 

3 If the density enhancement occurs in an angular-momentum- 
conserving fash ion, k oc S 1 / 2 and Qt oc S -1 / 2 as in the cas e of spi- 
ral shocks (e.g.,]BalbusXCowieI^ Rings 
in barred galaxies torm by gathering material that lost angular mo- 
mentum by a non-axisymmetirc bar torque as well as shocks. In 
our models, we found k calculated from the azimuthally-avcraged 
rotational velocity is similar to its initial profile, suggesting that 
density compression in the ring does not alter k much. 



tion. However, rings are already clumpy even in non-self- 
gravitating models as a result of the wiggle instability of 
dust-lane shocks. In models with intermediate bar torque 
(i.e, Qb ~ 0.1 — 0.2), therefore, clumps in nuclear rings 
simply acquire more mass and become denser due to self- 
gravity. In models with Qb XL 0.2, self-gravity produces 
additional fragments that make rings more clumpy in 
comparisons with non-self-gravitating counterparts (see, 
e.g., Figs. [2]{ and[3]<). The typical mass of clumps pro- 
duced is ~ 10 6 — 10 7 Mq. We note that the maximum 
density such clumps can attain is limited by numerical 
resolution. Self-gravitating clouds are r esolved only when 



Ar < Aj/4 (Truelove et al.|1997 1998), where Ar is the 



grid spacing and Xj = c z s /{GY>) is the local Jeans length. 
Therefore, clumps with £ ^ 10 3 M pc~ 2 at r ~ 1 kpc 
are gravitationally unresolved in our simulations. 

Another important effect of self-gravity is that it 
makes nuclear rings larger compared to those in non- 
self-gravitating models. Figures [8] and [10] show that 
compared to in Model M30R20N, the ring in Model 
M30R20G is larger in size at t = 0.3 Gyr and decays less 
afterward (see also Fig. S. This is because self-gravity 
deepens the total gravitational potential at the location 
of a ring, making the orbits of the ring material relatively 
intact and resistant to external perturbations that tend 
to change the orbits. Therefore, the radial decay of a ring 
due to the addition of low angular-momentum gas from 
outside is slower in self-gravitating models. The effect 
of thermal pressure in dispersing the ring material spa- 
tially is also smaller, resulting in a narrower and larger 
ring than in non-self-gravitating models. 

Finally, we discuss rings in self-gravitating models with 
a massive bar (/bar = 0.6). Figure [8] shows that Model 
M60R20G possesses a^-type, double rings at t = 0.3 Gyr, 
which is unlike Model M60R20N that has an inclined, sin- 
gle ring. As explained above, a large bar torque in Model 
M60R20N produces an inclined ring at t ~ 0.12 Gyr that 
processes with time to become of the xi-type. Strong 
self-gravity in Model M60R20G prevents the precession 
of the ring by providing additional non-axisymmetric 
torque, gradually aligning its long axis parallel to the bar 
minor axis. At t = 0.22 Gyr, the ring is still quite ec- 
centric and forms many high-density clumps via gravita- 
tional instability. When dense (unresolved) clumps move 
close to the galaxy center on their orbits (t ~ 0.26 Gyr), 
they gravitationally interact with other clumps at the op- 
posite side. These interacting clumps lower their orbits 
to smaller r, gather gas from the ring outside, and form 
a inner ring at t — 0.29 Gyr. The inner ring slowly dis- 
sipates as it loses mass through the inner boundary. In 
Model M60R25G, self-gravity results in an X2-type ring. 
But, the ring in this model is already so small that the 
gravitational interaction of dense clumps do not produce 
an obvious inner ring. Small, a;2-type rings in Models 
M60R30G and M60R35G (Fig.J%,t) become weaker with 
time and vanish at t ~ 0.75 Gyr, so that they are not 
considered as being permanent (Fig. (4Jd). 

4.3.3. Ring Properties 

To quantify the ring properties, we at each time calcu- 
late the mean density S r ; ng , mean radius r r ; ng , ellipticity 
firing = 1 — b r /a r with a r and b r denoting the semi-major 
and minor axes, respectively, and the position angle r i ng 
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Figure 12. Dependence on the bar strength Qi, of the ring size r r i ng relative (a) to the bar semimajor axis a and (b) to the ILR radius 
r lLRi ( c ) the ring ellipticity e r j ng , and (d) the ring position angle ipring f° r models with X2-type rings. In each panel, symbols give the mean 
values averaged over t = 0.3—1.0 Gyr, while errorbars represent the standard deviations. Open and filled sy mbo ls are for non-self-gravitating 
and self-gravitating models, respectively. Dotted and solid lines in (a) and (b) draw equations (j9| and |10| for non-self-gravitating and 
self-gravitating models, respectively. 



of the long axis of the ring, and take their temporal aver- 
ages over t = 0.3-1.0 Gyr. The resulting mean values and 
standard deviations are given in Tables 13] and H| for non- 
self-gravitating and self-gravitating models, respectively. 
Figure [12] plots these as various symbols and errorbars 
as functions of Qb for models with a^-type rings. Open 
symbols are for the non-self-gravitating models, while the 
self-gravitating results are plotted as filled symbols. It 
turns out that r r ; ng is best fitted solely by Qb not in 
combination with TZ. Our best fits of r,i ng against Qb 
are 



''rii 



0.062Q 
0.059Q 



-0.46 
-0.51 5 



without self-gravity, 
with self-gravity, 



(9) 



in terms of the bar major axis, and 



0.16Q b -°- 45 , 
0.19Q- 040 , 



without self-gravity, 
with self-gravity, 



(10) 



in terms of the ILR r adiu s. 

(|9l and ( 10 ) are plotted in Figure 12 i,b as 



Equations 

dotted and solid lines tor non-self-gravitating and self- 
gravitating models, respectively. Overall, rings in self- 
gravitating models are about ~ 5-20% larger than in 
non-self-gravitating models. The decrease of r r i ng with 
Qb is of course due to the fact that a stronger bar torque 
drives the gas closer to the galaxy center by removing 
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Table 3 

Properties of Nuclear Rings in Non-self-gravitating 



Models 



Table 4 

Properties of Nuclear Rings in Self-gravitating Models 



Model 



?/ s o 



(kpc) 



coring (deg) 



Model 



,/Eo r-ring (kpc) 



0ring (deg) 



M08R15N 
M08R20N 
M08R25N 
M08R30N 
M08R35N 



1.5 ± 0.2 
4.1 ± 1.5 
5.1 ± 1.8 

4.6 ± 2.0 



1.27 ± 0.19 
1.04 ± 0.21 
0.99 ± 0.19 
0.97 ± 0.14 



0.31 ± 0.05 
0.25 ± 0.05 
0.24 ± 0.05 
0.24 ± 0.04 



7.6 ± 4.7 
5.3 ± 4.3 
6.2 ± 4.4 
12.4 ± 3.4 



M15R15N 
M15R20N 
M15R25N 
M15R30N 
M15R35N 



1.5 ± 0.3 

3.6 ± 1.4 
9.3 ± 4.4 
9.8 ± 3.9 
14.1 ± 9.8 



1.33 ± 0.17 
0.97 ± 0.07 
0.81 ± 0.12 
0.76 ± 0.12 
0.72 ± 0.13 



0.32 ± 0.02 
0.20 ± 0.08 
0.19 ± 0.05 
0.18 ± 0.04 
0.20 ± 0.05 



8.3 ± 4.3 
8.8 ± 4.6 
17.2 ± 7.2 
31.9 ± 12.6 
36.9 ± 13.1 



M30R15N 
M30R20N 
M30R25N 
M30R30N 
M30R35N 
M60R15N 
M60R20N 
M60R25N 
M60R30N 
M60R35N 

Note. - 



2.5 EE 0.8 
8.2 ± 4.4 

22.9 ± 13.0 
28.9 ± 14.0 
38.0 ± 12.4 

4.6 ± 2.5 
1.4 ± 1.2 
0.8 ± 2.0 



-Ting j ' ring j 



1.00 EE 0.22 
0.74 ± 0.11 
0.60 ± 0.13 
0.56 ± 0.11 
0.42 ± 0.06 
0.63 ± 0.16 
0.15 ± 0.10 
0.21 ± 0.14 



and 



0.33 EE 0.06 
0.19 ± 0.06 
0.23 ± 0.06 
0.30 ± 0.14 
0.30 ± 0.10 
0.21 ± 0.09 
0.63 ± 0.22 
0.47 ± 0.17 



7.1 EE 5.6 
22.9 ± 9.3 

50.2 ± 17.2 
53.8 ± 23.2 

59.3 ± 24.5 
16.8 ± 13.5 
36.0 ± 36.7 
87.8 ± 17.3 



denote the mean density, 
mean radius, cllipticity, and position angle of the long axis of a ring, 
respectively. 

larger angular momentum from it. Note that r r j ng /riLR 
is less than unity and can be as small as ~ 0.2, demon- 
strating again that the ring position is not determined 
by the ILR but by Qb- We will show in Section [S] that 
the decreasing tendency of r r j ng with Qj, is qualitatively 
con sistent with the observational results of IComeron etl 
al. (2009) who found that nuclear rings are smaller in 



more strongly barred galaxies. In our models, nuclear 
rings are in general elliptical with e r j ng ~ 0.2-0.3, insen- 
sitive to Qb and self-gravity. They are well aligned with 
the bar minor axis when Qb < 0.1, while becoming more 
inclined as Qb increases. 

4.4. Nuclear Spirals 

There are two types of perturbations that can ex- 
cite nuclear spirals in the central regions: the non- 
axisymmetric bar potential and sonic perturbations from 
a nuclear ring. At early time before a ring forms, the 
bar potential is a lone perturbing agent of m — 2 spiral 
waves. As explained in Paper I, our galaxy models with 
a central BH have d(fl — K,/2)/dr > at r m j n < r < r max 
and d(Cl — re/2)/dr < otherwise, with the local mini- 
mum and maximum of the f2 — re/2 curve occurring at 
^min ~ 0.2 kpc and r max ~ 0.4-0.6 kpc depending on 
/bar and 1Z. Thus, the spiral waves at r < r m - m are 
trailing, while they are leading at r min < r < r max early 
time. Since the bar potential is nearly axisymmetric near 
the center, these waves would remain weak unless addi- 
tional perturbations are supplied. In models with an 
22-type ring, a ring that forms outside the spirals pro- 
vides sonic perturbations that propagate inward in the 
form of trailing waves. These waves interact construc- 
tively (destructively) with the trailing (leading) part of 
the spirals. Consequently, the outer leading parts are 
destroyed, while the inner trailing spirals grow stronger 
and extend outward to make contact with the ring (Fig. 
[8k,b). On the other hand, models with an Xi-type ring 
ao not retain nuclear spirals since the ring gas on eccen- 
tric orbits wipes out inner coherent spiral structures (Fig. 



M08R15G 
M08R20G 
M08R25G 
M08R30G 
M08R35G 



2.4 ± 1.0 
4.2 ± 1.5 
5.6 ± 1.9 
6.1 ± 2.2 



1.32 ± 0.23 
1.15 ± 0.18 
1.09 ± 0.15 
1.05 ± 0.15 



0.28 ± 0.02 
0.22 ± 0.05 
0.22 ± 0.04 
0.25 ± 0.03 



8.2 ± 5.1 
5.8 ± 3.3 
6.7 ± 4.3 
6.4 ± 3.2 



M15R15G 
M15R20G 
M15R25G 
M15R30G 
M15R35G 



1.6 ± 0.4 
5.6 ± 2.3 
9.2 ± 4.6 
12.6 ± 6.0 
111.5 ± 56.7 



1.39 EE 0.17 
1.00 ± 0.12 
0.93 ± 0.09 
0.85 ± 0.07 
0.75 ± 0.23 



0.32 EE 0.02 
0.17 ± 0.08 
0.19 ± 0.04 
0.20 ± 0.03 
0.29 ± 0.07 



8.3 ± 4.2 
6.6 ± 5.4 
11.3 ± 5.3 
12.3 ± 6.8 
14.6 ± 11.8 



M30R15G 
M30R20G 
M30R25G 
M30R30G 
M30R35G 
M60R15G 
M60R20G 
M60R25G 
M60R30G 
M60R35G 

Note. - 



2.8 EE 1.1 
11.2 ± 5.0 

18.4 ± 7.5 
81.2 ± 23.4 

78.7 ± 17.1 
5.2 ± 2.5 

32.8 ± 10.0 

69.5 ± 10.5 



1.04 EE 0.26 
0.86 ± 0.10 
0.71 ± 0.05 
0.60 ± 0.15 
0.50 ± 0.12 
0.67 EE 0.13 
0.47 ± 0.15 
0.35 ± 0.29 



and 



0.32 EE 0.07 
0.18 ± 0.07 
0.28 ± 0.03 
0.28 ± 0.11 
0.26 ± 0.13 
0.20 EE 0.08 
0.27 ± 0.07 
0.46 ± 0.20 



8.2 EE 5.1 
13.4 ± 9.3 
15.6 ± 5.4 
24.4 ± 16.5 
27.3 ± 23.0 
20.1 EE 15.1 
14.8 ± 8.4 
21.8 ± 25.3 



denote the mean density, 
mean radius, ellipticity, and position angle of the long axis of a ring, 
respectively. 

[8}:). Even though Model M30R35N has an x 2 -type ring, 
it is so small that strong perturbations from it inhibit 
the growth of spirals inside. Therefore, the presence of 
nuclear spirals requires two conditions: they should be 
of the X2 type and sufficiently large (perhaps larger than 

^max ) • 

Unlike dust lanes and nuclear rings, both of which re- 
main relatively stationary after the bar potential is fully 
turned on, we find that nuclear spirals do not achieve a 
quasi-steady state in that their peak density E pea k_and 
pitch angle i p vary substantially with time. Figure 13 i,b 



plots the temporal changes of X pca k and i p measured 
at r = 0.2 kpc for a few selected models. Dotted and 
solid lines correspond to non-self-gravitating and self- 
gravitating models, respectively. Compared to the other 
non-self-gravitating models, Model M30R25N forms nu- 
clear spirals earlier due to a larger bar torque. In gen- 
eral, nuclear spirals tend to unwind (i.e., i p increases) 
as they become stronger. This appea rs to be a generic 
conseq uence of the nonlinear effect. |Lee & Goodman 



( |1999[ ) showed analytically that dispersion relations of 
nonlinear waves involve the wavenumber, frequency, and 
wave amplitude all together, which is unlike in the linear- 
wave case where the wavenumber and frequency are in- 
dependent of the amplitude. They further showed that 
the wavenumber of nonlinear spiral waves is a decreas- 
ing function of the amplitude and that nonlinear trailing 
waves unwind and become more nonlinear as they prop- 
agate inward due to the increase in the angular momen- 
tum flux carried by the waves. This nonlinear growth 
and unwinding of traveling waves is consistent with the 
behavior of n uclear spirals formed in our models, 
also 



Figure 13 



shows that nuclear spirals in self- 
gravitating models unwind more slowly than the non-self- 
gravitating counterparts. This is presumably because 
non-axisymmetric waves corotating with the bar in self- 
gravitating models should have a larger radial wavenum- 
ber k r than those in non-self-gravitating models in order 
to have the same frequency. With larger k r , the an- 
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Table 5 

Time-averaged Mass Inflow Rate and Its Dispersion 



Model 



M 



AM 



Model 



M 



AM 





K lvi J 1 


-i-i 
1 


v ivi y 1 


) 




( vr 

V ivi y 1 


! 




) 


M08R15N 


2.4 x 10- 


3 


1.7 x 10- 


3 


M08R15G 


2.2 x 10" 


3 


1.6 x 10- 


3 


M08R20N 


2.3 x 10" 


3 


2.1 x 10- 


3 


M08R20G 


2.1 x 10" 


3 


1.9 x 10- 


3 


M08R25N 


2.9 x 10" 


3 


1.9 x 10- 


3 


M08R25G 


2.5 x 10" 


3 


1.7 x 10- 


3 


M08R30N 


3.1 x 10" 


3 


2.1 x 10- 


3 


M08R30G 


2.4 x 10" 


3 


1.5 x 10- 


3 


M08R35N 


3.2 x 10" 


3 


2.0 x 10- 


3 


M08R35G 


2.4 x 10- 


3 


1.3 x 10" 


3 


M15R15N 


2.3 x 10" 


y 


2.1 x 10" 


3 


M15R15G 


2.2 x 10- 


3 


2.0 x 10" 


3 


M15R20N 


3.2 x 10" 


3 


2.2 x 10- 


3 


M15R20G 


2.4 x 10" 


3 


1.8 x 10~ 


3 


M15R25N 


3.3 x 10" 


3 


2.3 x 10- 


3 


M15R25G 


3.0 x 10" 


3 


1.9 x 10- 


3 


M15R30N 


2.6 x 10 - 


3 


1.7 x 10" 


3 


M15R30G 


2.6 x 10" 


3 


2.3 x 10- 


3 


M15R35N 


2.5 x 10- 


3 


2.3 x 10- 


3 


M15R35G 


1.2 x 10- 


1 


7.7 x 10° 




M30R15N 


3.1 x 10- 


3 


2.7 x 10 - 


3 


M30R15G 


2.7 x 10- 


3 


2.2 x 10- 


3 


M30R20N 


5.2 x 10" 


3 


3.9 x 10- 


3 


M30R20G 


2.9 x 10- 


3 


2.3 x 10- 


3 


M30R25N 


2.0 x 10" 


3 


3.2 x 10- 


3 


M30R25G 


4.2 x 10- 


3 


2.9 x 10- 


3 


M30R30N 


2.5 x 10" 


2 


5.4 x 10- 


2 


M30R30G 


3.6 x 10- 


1 


2.0 x 10 1 




M30R35N 


1.7 x 10" 


1 


2.3 x 10- 


1 


M30R35G 


3.8 x 10" 


1 


1.9 x 10 1 




M60R15N 


7.5 x 10" 


3 


8.6 x 10" 


3 


M60R15G 


3.9 x 10" 


3 


4.6 x 10" 


3 


M60R20N 


3.2 x 10" 


1 


8.1 x 10- 


1 


M60R20G 


1.2 x 10- 


2 


2.6 x 10- 


1 


M60R25N 


5.4 x 10- 


1 


2.3 x 10° 




M60R25G 


6.3 x 10" 


1 


4.3 x 10 1 




M60R30N 


8.8 x 10" 


1 


2.2 x 10° 




M60R30G 


7.0 x 10" 


1 


3.2 x 10 1 




M60R35N 


9.2 x 10" 


1 


2.0 x 10° 




M60R35G 


1.2 x 10° 




5.5 x 10 1 




Note. - 


Time average 


of M is taken over t = 0.1 — 


1.0 Gyr. 
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Figure 13. Temporal evolution of (a) the peak surface density 
Speak an d (b) the pitch angle i p of nuclear spirals at r = 0.2 kpc in 
Models M08R15, M15R20, and M30R25. Dotted and solid lines are 
for non-self-gravitating and self-gravitating models, respectively, 
(c) Mean values (symbols) of i p averaged over t = 0.3-0.5 Gyr 
with the standard deviations (errorbars) as a function of Q;, for 
all models with nuclear spirals. Open and filled symbols are for 
non-self-gravitating and self-gravitating models, respectively. 



gle between the gas streamlines and nuclear spirals are 
smaller in self-gravitating models. This leads to a smaller 
departure of gas trajectories from the circular motions, 
resulting in lower £ pea k. 

All nuclear spirals in our models are logarithmic in 
shape. They grow with time and eventually develop 
into shocks, with the shock formation epoch delayed pro- 
gressively with decreasing Qt,. For instance, the spirals 
in Models M30R25N, M15R20N, and M08R15N become 
shocks at t ~ 0.26, 0.42, and 0.50 Gyr, respectively, after 
which Epeak decreases as most of the central gas inside 
the ring is accreted through the inner boundary. The 
shock formation time of nuclear spirals are delayed to 
t — 0.43, 0.45, and 0.55 Gyr for self-gravitating Models 
M30R25G, M15R20G, and M08R15G, respectively. At 
a given time, therefore, the nuclear spirals in smaller-Qh 
models tend to be more tightly wound, while those in 
larger-Qf, models tend to be more open and, sometimes, 
even shocked (Fig. [8]). This is illustrated quantitatively 
in Figure 13 z where the temporal averages of i p and the 
standard deviations over t = 0.3-0.5 Gyr are plotted for 
models with appreciable spirals during this time interval. 
The increasing behavior o f i n with Qb is consistent with 
the observational result of Peeples & Martini ( 2006 ) that 
nuclear spirals tend to be tightly wound in weakly barred 
galaxies. 

4.5. Mass Inflow Rate 

A galactic bar has often been considered to be a power- 
ful means to transport the interstellar gas at ~kpc scales 
all the way to the galaxy center, fueling AGN. To check 
the viability of this idea using our models, Figures 14 and 
\TE\ for non-self-gravitating and self-gravitating models, 
respectively, plot the mass inflow rates M through the 
inner boundary as functions of time. Table [5] gives the 
time-averaged values of M and the standard deviations 
AM over t = 0.1 — 1.0 Gyr for all models. One should be 
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Figure 14. Temporal variations of the mass inflow rates M for 
/bar = 0.3, and (d) / bar = 0.6. 

cautious in relating M to the accretion rate to a central 
BH, since the inflowing gas through the inner boundary 
may change its orbit due possibly to thermal and radi- 
ation pressures, gravity, etc. before reaching a BH and 
then come out of the inner boundary, or may be lost to 
star f ormation in a circumnucle ar disk surrounding a BH 
(e.g. " 




0.4 0.6 
t (Gyr) 

non-sclf-gravitating models with (a) /ij ar = 0.08, (b) /t, 



0.15, (c) 



Kawakatu & Wada 2008), which are not captured 



in our simulations. Therefore, M calculated from the 
current models can be considered as upper limits to the 
real accretion rates to a central BH. 

In non-self-gravitating models with a well-defined x-i- 
type ring, M is quite small at ~ 10~ 4 — 10 -3 M© yr -1 
at early time before the dust-lane shocks fully develop. 
It then increases as the gas flows inward along the dust 
lanes, but only to ~ 10~ 3 -10~ 2 M Q yr -1 unless the bar is 
strongly elongated with 1Z ^ 3.0. These relatively small 
values of M result from the fact that the majority of 
the inflowing gas from dust lanes is trapped in a nuclear 
ring that is located away from the center. Late time 
evolution of M depends on the size of a ring. In models 
with r r ing ~ 0.5 kpc, M remains below ~ 10 -2 M Q yr~ x . 
In Model M30R25N, on the other hand, the ring decays 
with time to have r r j ng ~ 0.4 kpc at t = 0.6 Gyr. At 
this small radius, thermal pressure is efficient to perturb 
a;2 _ hke gas orbits into eccentric shapes, increasing M to 



~ 1 M yr" 1 at t > 0.8 Gyr. 

For models whose central regions are dominated by 
rri-orbits, on the other hand, M increases dramatically 
to above ~ 1 M Q yr -1 , and sometimes as large as 
~ 10 M yr , as an inclined ring becomes smaller 
and more eccentric. When the short axis of an eccen- 
tric ring touches the inner boundary, the ring mate- 
rial flows in directly through it. This in turn causes 
the ring to decay rapidly in time, making M drop to 

- 10~ 2 - 10" 1 M yr" 1 . In Model M60R30N, gaseous 
blobs located in between the ring and the bar ends grad- 
ually move in toward the center and orbit along xi-orbits 
in the vicinity of the hole, which increases M again to 

- 1 M yr" 1 at t = 0.8 Gyr. 

In self-gravitating models with Qb 0.2, M is slightly 
smaller, owing to a larger nuclear ring, than that in 
the non-self-gravitating counterpart. In models with 
Qb £ 0.2, however, self-gravity makes the rings unstable, 
producing small dense clumps. While these clumps move 
along the rings, they interact with each other and some- 
times plunge into the central hole directly, increasing M 
instantaneously. In Model M15R35G with Q b = 0.21, 
such events occurring three times over 1 Gyr lead to 
M > 10 2 M yr -1 . In other models with larger Qi, such 
as Models M30R30G and M60R35G, direct accretion of 
dense clumps occur much more frequently, making the 
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Figure 15. Same as Figure fl4] but for self-gravitating models. 

M curves vs. time highly intermittent. 

5. SUMMARY AND DISCUSSION 

We have presented the results of high-resolution hy- 
drodynamic simulations on the formation and evolution 
of gaseous substructures in barred galaxies with varying 
bar strength. We initially consider an infinitesimally- 
thin, isothermal, unmagnetized gas disk with uniform 
surface density embedded in the external gravitational 
potential. We run both non-self-gravitating and self- 
gravitating models, but the effects of star formation and 
feedback are not considered in the present work. In or- 
der to focus on the effects of the bar parameters, we fix 
the gas sound speed to c s = 10 km s _1 and the BH 
mass to Mbh = 4 x 10 7 M Q that affects the rotation 
curve near the galaxy center, and vary two parameters: 
the bar mass measured by /b ar (Eq. pi) relative to the 
spheroidal component and its aspect ratio TZ (see Tables 
[l] and [2] for model parameters) . 

In what follows, we summarize the main results of the 
present work and discuss them in comparison with ob- 
servations. 

1. Bar Strength Parameter - We measure the bar 
strength using the dimensionless parameter Qb (Eq. [4]) 
as is usually done in observational studies. For our galaxy 
models with a Ferrers prolate bar with index n, we calcu- 
late Qb as well as the radius r{Qb) where the maximum 
bar torque occurs, and provide the fitting formulae (Eqs. 



u 

© 

a 
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|5] and [6]) as functions of /bar and TZ. While Qb is 
linearly proportional to TZ and almost linearly to /ban 
r(Qb) is a weakly decreasing function of 1Z and indepen- 
dent of /bar- Both Qb and r(Qb) become smaller for a 
more centrally-concentrated bar. 

Having found the dependence of Qb and r{Qb) on the 
other bar parameters in our galaxy models, it is inter- 
esting to apply our results to o bserved barred galaxies. 
Recently, Comeron et al. (20101 measured Qb, r(Qb), TZ 
(or, equivalently, the bar ellipticity) , and the bar semi- 
major axis a for a sample of nearby galaxies that contain 
nuclear rings, which is by far the most complete sample. 
Figure [16k plots as sta r symbols the empiric al relation 
betweenQb and TZ from Comeron et al. (2010) for galax- 
ies with 1.5 < 72. < 3.5^ Also plotted are equation ([5| 
for various values of n and /bar- Note that the trend of 
Qb becoming larger for larger TZ in the observational esti- 
mates is entirely consistent with the results of our galaxy 
models. 

When approximating the observed bars using Ferrers 
prolate spheroids, the observational results for Qb vs. 
TZ can be best described by the bar mass fraction of 
/bar = 0.3-0.5 for n = 1 inhomogeneous bars and 
/bar = 0.25-0.35 for n — homogeneous bars. The ob- 
served relation between r(Qb)/a and TZ shown in Fig- 
ure |16b appears better explained by the n — homo- 
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Figure 16. Relationships (a) between Qi, and TZ = a/6 and 
(b) between r(Q(,) and TZ from our galaxy models (various lines) 
with a Ferrers prolate bar (Eqs. [5] and [6]) in comparison with 
observational results (star symbols; oflComeron et al.l (12010b . The 
observed bars are best represented by /bar = U.t> and ri < 1. 

geneous bars than inhomogeneous barsj^] We note that 
many uncertainties surround the observational determi- 
nations of Qb and r(Qb) as they rely sensitively on the 
bulge subtraction, assumptions on the disk scale height 
and orientation angle, etc., which are quite uncertain 



(e.g., Laurikainen et al. 



Given that observationa 



2004[ |2006| |Buta et al |2006| ) 

errors involved m the Qb (also 
likely in r(Qb)\ deter minations are typically ~ 20% (e.g 



Comeron et al.|2009[ ), the comparison shown in Figure 16 
suggests that bars in real galaxies are most likely to have 
a mass fraction /bar = 0.25-0.5 of the spheroidal com- 
ponent, and unlikely to be more centrally concentrated 
than the n = 1 case. 

2. Dust Lanes - The imposed non-axisymmetric bar 
potential readily induces dust-lane shocks across which 
gas in rotation about the galaxy center loses angular 
momentum significantly and falls radially inward. Dust 
lanes in a quasi-steady state approximately follow one of 
a^i-orbits aligned parallel to the bar major axis. The cur- 
vature Aa of dust lanes in our models depends primarily 
on Qb in such a way that they tend to be more curved 
under a weaker bar. This results from the facts that dust 
lanes are closer to the bar major axis as Qb increases and 
that inner xi-orbits are more eccentric than outer ones 

4 Galaxies with r(Qb)/a > 1 have non-axisymm etric torques 
dominated by outer spiral arms rather than by bars | |Comer6n et| 
|al.|2010| . 



Figure 17. Comparison of our numerical results (squares with 
errorbars) with the observational measurements for the relations 
(a) between the dust-lane curvature Act and Qt, and (b) between 
the ring radius r r ; ng and Qb- Open and filled squares are for non- 
self-gravitating and self-gravit ating models, respectiy ely. In (a), 
star symbols are adopted fromlComeron e t al,| ||2009|, while star 
sym bols and filled triangle s in (bj are trom (Jomcron et al.| l |2010| l 
and |Mazzuca et al,| ( |201l} , respectively. 

(see Fig. [7]). It also depends, albeit less sensitively, on 
1Z in that a more elongated bar has more straight dust 
lanes (Eq. [8]), since xi-orbits are rounder with smaller 
7Z when Qb is fixed. Dust lanes are not much affected by 
self-gravity since they have a Toomre stability parame- 
ter greater than unity and are characterized by strong 
velocity shear. 

As mentioned earlier, that dust lanes are more straight 
under a s tronger bar potential w as first theoretically pre- 
dicted by Athanassoula (1992b I and lat er confirmed em- 
pirical ly by" Knapen et al. ( 2002 ) and iComeron et al 



(120091) . In particular, Comeron et al. 
the dust-lane curvatures in a sample o 



measured 
barred galax- 



ies that contain clear dust lanes in the SDSS DR7 or 
NED images, and studied a relationship between Aa and 
Qb- Figure [17^ reproduces their observational results as 
star symbols in comparison with our numerical results 
shown as open and filled squares with errorbars for non- 
self-gravitating and self-gravitating models, respectively 
Note that Aa decreases with increasing Qb in both nu- 
merical and observational results, although Aa in the 
simulations corresponds roughly to the lower envelope of 
the observational results. Note also that our numerical 
results are unable to reproduce the scatter seen in the 
observations. These quantitative differences of Aa be- 
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twe en ou r numerical work and the work of IComeron etl 
al. (2009) are likely due to the differences in t he spatial 
range s or dust lanes where Aa is measured. IComeron 



rings smaller (Pa per I) 



et al. ( 2009 1 considered a constant-curvature range that 
varies from galaxy to galaxy, while we fix the position 
angles of both ends of the range. Also, the fact that 
the parameter space covered by our models is very lim- 
ited (i.e., fixed sound speed, n = 1 Ferrers bar potential, 
no magnetic field, etc.) may also be partly responsible 
for the differences between our numerical results and the 
observations. 



More recently, 



Comeron et al. 



(2009) further noted that there is a 
large spread of Aa tor given Qb- To find the origin of 
the spread, they ran a large number of numerical simula- 
tions with differing bar and bulge parameters, and found 
that the spread in Aa can be reduced if the aspect ra- 
tio TZ is considered together with Qb when fitting the 
curvatures. We similarly found that a linear combina- 
tion of Qb and TZ provides a better fit than Qb alone, 
although Aa decreases with increasing TZ in our fit (Eq. 
[8]), while Aa is an increasing fu nction of TZ in their fit 



[3] of Comeron et al. 2009). This discrepancy is 



again thought ot as arising from the differences in the 
ranges of dust lanes where Aa is measured and from the 
limited range of the parameter space in our models, as 
mentioned above. 

3. Nuclear Rings: Size - The shocked gas moving in 
toward the central regions along the dust lanes has so 
large a speed that the bar torque cannot stop its mo- 
tion across the ILR. The inflowing gas keeps moving in 
and eventually forms a nuclear ring at the location where 
the centrifugal force balances the external gravitational 
force. The mean radius of a ring in our models is gener- 
ally smaller than the ILR location, and decreases system- 
atically with increasing Qb (Eqs. [9] and 10 ). By making 
the total gravitatioal potential deeper, self-gravity makes 
the ring larger in size by ~ 5 — 20%, with larger values 
corresponding to lager Qb- Combined with the results of 
Paper I that showed that the ring position is insensitive 
to the mass of a central BH (and thus to the number 
and locations of ILRs), this clearly evidences that the 
ring position is not determined by the resonant interac- 
tions of gas with the underlying gravitational potential 
but rather by the amount of angular momentum loss at 
the dust-lane shocks. 



Comeron et al. (20101 also presented the sizes of nu- 
clear r ings in their sample of barred spiral galaxies. 



Figure 17 d compares their results (star symbols) with 
our simulation outcomes (squares with errorbars) on the 
fring/a-Qb plane. Both observational and numerical re- 
sults show that r r ; ng /a becomes smaller with increasing 
Qb, indicating that stronger bars can possess smaller 
rings. For Qb 0.15, the agreement between observa- 
tional and numerical results is quite good. For Qb Si 0.15, 
on the other hand, r r i ng /a in our models corresponds 
roughly to the upper envelope of the observational re- 
sults for given Qb- This is presumably because our nu- 
merical models are unmagnetized: inclusion of magnetic 
fields efficiently removes angular momentum further at 
dust-lane shocks, which makes the ring size smaller by a 
factor of ~ 2 when magnetic fields have an equipartition 
strength with the thermal energy (Kim & Stone 2012). 
A larger effective sound speed can additionally make the 



(2011) measured the 



Mazzuca et al 
ring radii for a sample of 13 barred/ unbarred galax- 
ies that contain star-forming nuclear rings: their results 
for 8 gala xies whose bar sizes a re g iven in | Comer on et 
al. J 2010 ) a re plo tted in Figure |17b as triangles. |Maz- 



zuca et al. (2011) argued that the ring size is well cor- 



1S 



related with the compactness C = v$/r t , where r t 
the turnover radius of the rotation curve that has the 
velocity vq at the fiat part, such that more compact 
(with smaller turnover radius) galaxies have a smaller 
ring. In our galaxy models, however, r t ~ 1 kpc and 
C ~ 4 x 10 4 (km s -1 ) 2 kpc -1 , insensitive to Qb and TZ, 
while the variation of the ring size with Qb is about by 
a factor of ~ 2-3. Since the sample galaxies in [Mazzuca 



et al. (2011 ) also exhibit a positive correlation between C 
and Qb, the negative correlation between C and r r i ng in 
their results may simply be a reflection of a more intrin- 
sic negative correlation between r r ; ng and Qb- Of course, 
our current models with a limited range of the parame- 
ters cannot address the effects of C on the ring size. 

4- Nuclear Rings: Shape - Not all the rings have a con- 
ventional shape similar to X2-orbits in our models unless 
self-gravity is included. In non-self-gravitating models 
with /bar = 0.6 and TZ > 2.0, the inflowing gas can- 
not settle on a^-orbits since they have too small kinetic 
energy. The ring material instead takes on an inclined 
orbit in between x\- and xi -families. An inclined ring 
in these models becomes smaller and more eccentric as 
low angular-momentum gas is added from outside. It 
loses much of its mass through the inner boundary when 
its short axis moves close to the center. Due to the bar 
torque, the inclined ring precesses slowly to align its long 
axis parallel to the bar major axis, eventually forming an 
xi-type ring. Strong self-gravity of an inclined ring pro- 
vides additional non-axisymmetric torque that prevents 
the precession of an inclined ring, leading instead to an 
^2-type ring in self-gravitating models. 

When /bar < 0.3 or TZ < 1.5, on the other hand, the 
a^-family of closed orbits near the center have sufficiently 
large kinetic energy that the inflowing gas along the dust 
lanes can transit easily to one of them in both non-self- 
gravitating and self-gravitating models. These rings are 
eccentric with an ellipticity of e r i ng ~ 0.2-0.3 insensitive 
to Qb, which is within the rang e of the observed ring 
cllipti cities, ~ 0-0.4 , repor ted by Comeron et al. (2010) 



and Mazzuca et al. (2011). Again, magnetic fiel ds ar e 
expected to circularize nuclear rings (Kim & Stone 2012). 

5. Nuclear Spirals - Well-defined twin-armed nuclear 
spirals grow only in models in which the nuclear ring is 
of the X2-type and sufficiently large in size: they would 
otherwise be destroyed by the ring material on highly 
eccentric orbits. Even in models with an X2-type ring, 
nuclear spirals are absent if the ring is so small to limit 
their spatial extent. All nuclear spirals in our mod- 
els are trailing and logarithmic in shape. While the 
shapes of dust lanes and nuclear rings do not change 
much with time after the potential is fully turned on, 
nuclear spirals are not stationary over the course of the 
entire evolution. They initially start out as being tightly 
wound and weak, and then gradually unwind and become 
stronger until turning into shocks. This unwinding and 
growth of nuclear spirals appears to be a generic prop- 
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erty of nonlinear waves that become more no nlinear as 
they propagate inward ( Lee fc Goodman|1999 ) . The un- 
winding rate is lower in self-gravitating models than in 
non-self-gravitating models. Since nuclear spirals grow 
and unwind faster as Qb increases, the probability of 
having more tightly-wound and weaker spirals is larger 
for galaxies with a weaker bar torq ue. This is consis- 
tent w ith the observational results of |Peeples fe Martini 
(2006) who found that tightly wound spirals are found 



primarily in weakly barred galaxies, while loosely wound 
spir als are more common in st rongly barred galax i es (see 
also |Martini et al.|[2003a|b[ ). |Peeples fc Martini| p006 | 
also found that grand-design nuclear spirals in strongly 
barred galaxies does not extend all the way into the nu- 
cleus, which is consistent with our numerical results that 
show that nuclear spirals cease to exist by turning into 
shocks, and this happens earlier in highcr-Qf, models. 

6. Mass Inflow Rates - In our models, the mass in- 
flow rate to through the inner boundary is found to be 
M ~ 10~ 3 - 10~ 2 M yr" 1 for models with Q b < 0.2 
regardless of the presence of self-gravity. Without self- 
gravity, models with Qf, ^ 0.2 but with an X2-type ring 
still have M < 10 -2 M Q yr _1 , since most of the inflow- 
ing gas is trapped in the ring. If the inflowing gas moves 
all the way to a central BH, these values of M corre- 
spond to the Eddington ratio A = iboi/^Edd = 4.5 x 
10- 2 (Af/10- 2 Mq ^"^(Mbh/IO 7 Mo)- 1 - lO^lO" 2 
(Paper I), pote ntially exp laining low-luminosity Seyfert 
1 AGNs. (e.g., |Ho||2008 [). Here, Lboi and ^Edd denote 
the bolometric and Eddington luminosities of an AGN, 
respectively, and 10% of the mass-to-energy conversion 
efficiency of the accreted material is assumed. Some of 
our numerical models exhibit unrealistically large mass 
inflow rates. Non-self-gravitating models in which the 
central regions are dominated by the gas on x\ -orbits 
are found to have M ^ 1 M Q yr _1 and sometimes as 
large as ~ 10 M Q yr" 1 when the ring gas on eccen- 
tric orbits is accreted directly to the inner boundary. In 
self-gravitating models with Qi, ^ 0.2, on the other hand, 
rings are unstable to form high-density clumps with mass 
~ 10 6 — 10 7 M . These clumps sometimes plunge into 
the central hole, causing the mass inflow rate to fluctuate 



with large amplitudes. 

We finally remark some caveats associated with M 
obtained in our simulations and in interpreting it as 
a mass accretion rate to a central BH. First, M from 
sclf-gravitating models are definitely more realistic that 
that from non-selfOgravitating models. Still, Models 
M60R30G and M60R35G can not be applied to real 
galaxies since they do not posses well-defined nuclear 
rings: the bar in these models is perhaps too massive 
or too elongated, or the Ferrers prolate bar is not a 
good representation of realistic bars. Second, as men- 
tioned earlier, M measured is the rate of gas mass that 
goes in through the inner boundary. This is likely to 
be an upper limit to the real accretion rate to the BH 
since some of the inflowing mass changes its orbit before 
reaching the BH and can possibly come out of the in- 
ner boundary. Third, high-density clumps produced by 
gravitational instability of nuclear rings would undergo 
star formation, reducing gas content in the rings. Ensu- 
ing feedback would destroy them, so that the accretion of 
dense clumps in real situations would much less frequent 
than in our simulations. Fourth, a circumnuclear disk 
with starburst activities can make t he accretion rate to 
the BH much smaller than M (e.g, Davie s et al. 2007 
Watabe et aT1|2008| |Kawakatu fc Wada| |2008[ ) . Fifth, 



while magnetic fields are known to enhance M consider- 
ably ( |Kim fc Stone|2012[ ), they would suppress or reduce 
gravitational instability of nuclear rings, which tends to 
reduce M. In order to properly evaluate the mass accre- 
tion rates, therefore, it is required to run more realistic 
models of barred galaxies including star formation, feed- 
back, magnetic fields, and other physical processes that 
affect gaseous features and gas inflows in the nuclear re- 
gions. 
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APPENDIX 

BAR STRENGTH FOR A PROLATE FERRERS BAR 

In this Appendix, we provide analytic expressions for the bar strength Qt and the radius r(Qb) of the maximum 
bar torque for a galaxy with an n = 1 prolate Ferrer bar, whose density distribution is given by equation Q. The 
eccentricity of the bar is e = (1 — 6 2 /a 2 ) 1//2 . Since the maximum bar torque always occurs inside the bar, it is sufficient 
to consider the interior (i.e., g < 1) gravi tational potentia l of the bar. 



Following the procedures presented by Pfenniger (1984) for triaxial ellipsoids and by Maciejewski & Sparke (20001 
for prolate spheroids, one can show that the interior potential of a Ferrers bar with n 
by 



1 at the z = plane is given 
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(A6) 

(A7) 



= -2irGab 2 p ha , r xy{C x x 2 + C + C y y 2 )/r, 



Note that /, A\, and ^3 in equations (A2) - (A4) are i dentical to the dimensionle ss coefficients for the gravitational 
potential of uniform spheroids tabulated in Table 2.1 of Binney & Tremaine (2008). 
The gravitational force in the azimuthal direction due to the n = 1 prolate is then given by 

1 9$ bar 

±<T = 

r c 

where 

C x =W 02 -W n , 
C = W 10 -W 01 , 
C y = W n - W 20 . 



(A8) 



(A9) 
(A10) 
(All) 



It can be easily verified that C x , C y > and Cq < for < e < 1. In the limit of e — > 0, C x = C y — > 2e 2 /(7a 5 ), while 
C -> -2e 2 /(5a 3 ). 
It is straightforward to show that the maximum value of \Ft\ is 



47rGa6 2 p bar (-C0) 3 / 2 



3 3 / 2 



1/2 



c. 



1/2 ' 



(A12) 



which occurs at x 2 = -C /(3Cy 2 )(Cy 2 + C^ 2 )- 1 and y 2 = -C /(3C y /2 )(C x /2 + Cy 72 )" 1 . 

For galaxies with a constant rotational velocity Vq, the radial gravitational force is Fr = v 2 /r. Therefore, the bar 
strength parameter defined in equation Q becomes 



which is attained at x 



Q b = 

-C /(AC X ) and y 2 



(rF T ) r 



7rGa6 2 / 9 b ar Cl 



Av 2 



^JC X Cy 



-Co/(4C a ), or at the radius 



r(Qb) 



(-Co) 1 / 2 / 1 



C x Co 



1/2 



(A13) 



(A14) 



Let Mtot(a) = clvq/G denote the total galaxy mass within r < a. Assuming that the galaxy mass is dominated by the 
bulge and the bar inside r = a, it then follows that 



_ 15 aC 2 



(A15) 



Figure [I] plots Qb and r(Qb) in equations (A15) and (A14| as dotted lines, which is not much different from the true 
bar strength that does not assume flat rotation. 
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